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We construct and analyze climate networks based on daily satellite measurements of temperatures 
and geopotential heights. We show that these networks are stable during time and are similar over 
different altitudes. Each link in our network is stable with typical 15% variability. The entire 
hierarchy of links is about 80% consistent during time. We show that about half of this stability is 
due to the spatial 2D embedding of the network, and half is due to physical coupling mechanisms. 
The network stability of equatorial regions is found to be lower compared to the stability of a typical 
network in non-equatorial regions. 
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I. INTRODUCTION 

During the past decade, methods from network the- 
ory have been applied to describe complex systems that 
are composed of many interacting components (see e.g. 
[l|-[26|). While in some cases the representation of such 
systems as a network is obvious and the nodes and links 
are identified directly (e.g. cables connecting computers 
in a computer network) |2| , in many real world networks 
such as biological systems [1], neural networks [1^, cli- 
mate networks |27H35| and others, the identification of 
links is not direct and it is based on statistical analysis 
of the similarity of the dynamics of nodes [36j . 

In recent years it was suggested that climate fields such 
as temperature and geopotential height at a certain pres- 
sure level can be represented as a climate network where 
the nodes are geographical sites and the links are the 
information fiow between these sites (nodes). Although 
the dynamics within a single node is unpredictable and 
chaotic, yet the dynamic of one node may be coupled to 
the dynamics of other nodes and could be observed. The 
correlations between the dynamics in two different nodes 
is represented in our network as a link between them. 
Since we deal with climatological data, the correlations 
(Sec. Ill Bp might be with a time-delay. Therefore, each 
link is quantified by two parameters; a strength value 
which quantifies the intensity of the cross correlations, 
and a time delay value which quantifies the delay in the 
data flow between the two nodes. 



II. THE CLIMATE NETWORK 
A. Data 

We analyze data obtained from a reanalysis 
project [131 • The records consist of the reanalysis 
air temperature field and the geopotential height field. 
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for the lOOOhPa, 925hPa, 850hPa, TOOhPa, 500hPa and 
SOOhPa isobars. We use daily values between the years 
1948-2006. The data is arranged on a world-wide grid 
with a resolution of 5° x 5°. We divide the globe into 9 
zones (see Fig. [1]), in order to identify different network 
dynamics specific to different zones. 




FIG. 1: (color online). The geographical locations of the 9 
separate zones, on which we base our network analysis. 



B. The network construction method 

We analyze daily climatological records (tempera- 
ture/pressure) taken from a grid in various geographi- 
cal zones (Fig. [1]). To avoid the trivial effect of seasonal 
trends we subtract from the records of each day the yearly 
average value of that day. Specifically, we take the cli- 
matological signal (temperature/pressure) of a given site 
in the grid to be S^{d), where y is the year and d is the 
day (ranging from 1 to 365) of that year. The new sig- 
nal will be Sy{d) = Sy{d) - ^ Ey S^id), where N is the 
number of years available in the record. For each pair 
of sites I and r in a specific zone, we compute the abso- 
lute value of the cross-covariance function Xf^ of their 
local climatological signals such as temperature/pressure 
in the range of time delays r G [—Tmax, Tmax] integrated 
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link to be 

where < 



over a specific year (y) (see Fig. [2|). Using the cross- 
covariance function X^^ we define the strength of the 
{MAX{Xl)- < XI >)/STDiXl), 
>, MAX and STD are the mean value, 
maximal value and the standard deviation of Xj^^ in the 
range of r, respectively. The matrix Wj^^ represents the 
weighted adjacency matrix of the network at year y. The 
time shift at which Xf^ is maximal is defined as the link 
time delay, and denoted as T^'^. 




FIG. 2: A typical cross-covariance function between two sites, 
representing the level of correlation within a time lag ranging 
from -Tmax to +Tmax where Tmax = 220. In the current 
example WJ"^ = "■=17°,°^ = 5, Tj"^ = -15. 



III. RESULTS 



A. Stability of single links 



with time [38l - l43j . Therefore, this stability of the links 
over many years is surprising and may suggest that we 
can extract new information from the links between the 
nodes. 
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First we focus on the behavior of single links. We an- 
alyze the yearly variations of the link strength Wj^^ . We 
find that a typical link maintains its strength Wi_r during 
the years with typically small fluctuations of 15%. This 
stability in the strength of the links behavior is valid for 
links across long and short distances (see Fig. [3]). In 
Fig. 13] we show the distribution of the relative standard 
deviation of all links in two networks (a) the network lo- 
cated at Zone 1 (a non-equatorial region) and (b) the 
network located at Zone 9 (an equatorial region) (see 
Fig. [T]). We observe a significant difference in this dis- 
tribution between networks located in equatorial regions 
and networks located in non-equatorial regions. While 
in networks located in equatorial regions (zones 7-9), the 
minimum variation is about 0.1 and the maximum is 
about 0.3, in networks located in non-equatorial regions 
(zones 1-6), the minimum variation is about 0.05 and 
the maximum is about 0.25 (see Fig. |3|). Therefore the 
link strengths, Wi^r, in non-equatorial regions tend to be 
more stable than the link strengths, Wi^r, in equatorial 
regions. 

A typical auto-correlation function for climatological 
records of a specific node decays rapidly as a power low 



FIG. 3: Two examples of typical dynamics of a link strength 
(W) during the years. In (a) the distance between the two 
sites is about 750Km, the average time delay T — (day) and 
the variation in the link strength, STD{Wi^r)/Wi^r = 0.1. In 
(b) the distance between the two sites is about 1500-K'm, the 
average time delay T = 1 (day) and the variation in the link 
strength, STD{Wi,r) = 0.1. 



Analyzing the influence of spatial distances D, between 
the nodes on the strength Wd of the link leads to the ob- 
servation of a strong dependence of Wd on D (Fig. [5|). 
Here Wd = W^,^ is the average over all link strengths W 
at distance D, and over all years, y. It is seen that for 
D > 2000 Km, Wd reaches a low and almost constant 
value. This constant value can be regarded (as will be 
seen in Sec IIV Al) as the level of noise. However, we ob- 
serve a significant difference in this dependence between 
networks located in equatorial regions and networks lo- 
cated in non-equatorial regions (compare Fig. 
Fig. 



5(a) and 



5(b)). In networks located in equatorial regions 



(zones 7-9), Wd decreases significantly slower with time 
compared to other regions. This difference mainly ap- 
pears in networks based on the geopotential height field. 
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FIG. 4: Distribution of the links variation during time, 
STD (W) /W in two zones. Zone 1 and Zone 9, botii for 
network based on temperature at 850hPa isobar. 



B. Stability of the entire network 



In Sec. nil Al we showed that single Unks remain rela- 
tively stable. In this section we study the stability of the 
entire hierarchy of links within the climate network. We 
find the network to be relatively stable over time. This 
stability of the network may be demonstrated by mea- 
suring the similarity between network states in different 
years. We analyze the similarity by calculating pir) (the 
Pearson coefficient) between the adjacency matrices of 
two network states in different years yi ,2/2 as a function 
oi T = y2 — Ui- In Fig. [6] (the upper curve) we show the 
average similarity p(r) between network structures as a 
function of the time separation r. It is seen in Fig. 6(a 



that this similarity is indeed high and almost constant, 
p{t) 0.8. This behavior is consistent for all networks 
in the non-equatorial regions. For networks in equatorial 
regions the correlation between the network states in dif- 
ferent years, p{t) is still significantly high but relatively 
fluctuative. (see Fig. |6(b)| . 

In Fig. [5] we showed, that there is a strong dependence 
between the link strength, W and the link distance, D. 
Links with shorter distances D are therefore more likely 
to have higher W values, at all times. It is therefore plau- 
sible that the high stability of |5(t) « 0.8 is partially due 
to this strong dependence. The contribution of the effect 
of the W-D dependence to this observed stability, on 
the one hand, and the contribution of physical coupling 
processes, on the other hand, must be estimated. 

We achieve this goal of removing the contribution 
of the W-D dependence by subtracting from each link 
strength, Wi^r the average strength of the group of links 
with a similar distance, Wd- A new, transformed adja- 
cency matrix, Wi^r — Wd, is thus formed, which does not 
depend on D. Repeating our analysis of calculating p(r), 
for the new adjacency matrix, we show in the lower curve 
of Fig. [6l the stability of this network. Indeed, after re- 
moval of the distance effect, the network exhibits lower 
p(t) values. However, the stability related to physical 
coupling processes is still significant [47]. 
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FIG. 5: The dependence of Wi^r on Di.r in two typical loca- 
tions, (a) Zone 1 (non-equatorial region) . (b) Zone 9 (equato- 
rial region). The four curves describe four networks which are 
based on geopotential height measurement at 850 and 500hPa 
isobar and based on temperature measurement at 850 and 
500hPa isobar. 



It is plausible that some of our network links emerge 
mainly due to noise and not due to real physical coupling 
processes. In Sec. IIV Al we show that this group of false 
links is characterized by low Wi_r values at all times, and 
that this characterization is sufficient for uniquely iden- 
tifying this group. In order to identify this group of false 
links which are due to noise, we define Wi,r as the aver- 
age strength of a link Wl'^. over all years. Upon eliminat- 
ing low weighted links that satisfy Wi^r < & (where 6* is a 
threshold that will be determined later in Sec. lIV A| from 
our network, we observe in the middle curve of Fig. [6] an 
increase of the network stability. Thus, while the hierar- 
chy of false links rapidly changes in each time step, the 
hierarchy of significant links (having Wi^r > 0) is, to a 
large extent, preserved. 

The new p{t) values, which are calculated after the re- 
moval of both the distance effect and the effects of noise, 
for different regions and fields are summarized in Table 
|T1 In contrast to the common p{t) = 0.8 value that was 
observed for the original network (including the distance 
and noise effects), after removal of the distance and noise 
effects we observe lower p{t) values specific to each zone 
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FIG. 6: The average correlation, p(r), between network ad- 
jacency matrices at different time snapshot, between r = 1 
and r = 40 years apart, for networks based on temperature 
at 850hPa and located at (a) zone 1, and (b) zone 9. The 
upper curve in each figure, represents the correlation between 
the original networks without removing the effects of distance 
and noise. The lower curve represents pij) for a network af- 
ter removing of the distance effect. The mid curve represents 
P(t) for a network after removing both, distance and noise 
effects. 



and climate variable. Still, in general, non equatorial 
regions exhibit larger stability values than equatorial re- 
gions. 

It has been shown that during El-Niho times, link 
strengths, Wj'^, are significantly reduced mainly in equa- 
torial regions [13: [IE]- Hence our observation of lower sta- 
bility in equatorial regions is consistent with the known 
effect of El-Niho on the climate network. 

From Table U we see that removing both the distance 
effect and the effects of noise reveal that the networks in 
zone 3, at the southern ocean, exhibit low stability, p{t) 
values similar to the equatorial regions. This similarity 
of the behavior of the network in zone 3 and the behav- 
ior of the network in equatorial regions is consistent with 
the known local oscillations in zone 3 that correlate with 
ENSO, due to both ocean mechanisms [44] and atmo- 
spheric coupling mechanism j45l] . 



zone 

climate^^ 
variable ^^^^ 


Temp 
850 hPa 


Temp 
500 hFa 


Height 
ooO hra 


Height 
500 hra 


zone 1 


0.42 ± 0.08 


0.31 ±0.08 


0.48 ± 0.08 


0.3 ±0.1 


zone 2 


0.4 ±0.08 


0.33 ± 0.08 


0.32 ±0.11 


0.3 ±0.11 


zone 3 


0.27 ±0.09 


0.21 ±0.1 


0.17 ±0.13 


0.18 ±0.11 


zone 4 


0.32 ±0.08 


0.24 ±0.09 


0.45 ± 0.11 


0.36 ±0.1 


zone 5 


0.42 ± 0.1 


0.34 ± 0.1 


0.51 ± 0.12 


0.42 ± 0.11 


zone 6 


0.43 ±0.08 


0.41 ±0.09 


0.45 ± 0.08 


0.37 ±0.08 


zone 7 


0.22 ±0.11 


0.37 ±0.15 


0.37 ±0.16 


0.18 ±0.17 


zone 8 


0.33 ±0.13 


0.19 ±0.14 


0.24 ±0.15 


0.14 ±0.21 


zone 9 


0.25 ±0.14 


0.2 ±0.11 


0.37 ±0.14 


0.16 ±0.2 



TABLE I: The average correlation values, |5(r), between net- 
work adjacency matrices at different time snapshot for net- 
works based on various fields and located at different regions. 
The values shown are after removing the distance and noise 
effects. 



Based on the high stability values seen in Table HI 
we conclude that similarity between network states at all 
times stems from a hierarchy of real physical correlations 
(links) between different locations, which is preserved in 
time. 



IV. SIMILARITY BETWEEN THE NETWORKS 
STRUCTURE, IN DIFFERENT ALTITUDES AND 
DIFFERENT CLIMATE VARIABLES 

A further indication that the stability of the network 
structure reflects a stability of physical coupling pro- 
cesses, is from the finding of similarity between the net- 
works structure in different altitudes, and different cli- 
mate variables. For example, synchronized heating of 
two sites at the 850hPa isobar network is likely to also 
cause synchronized heating of the corresponding sites in 
the adjacent isobar of 500hPa network by direct heat 
transport. In the following we show that such a corre- 
spondence between networks of different altitudes and 
climate variables indeed exists. 

In Fig. [7] we show the Pearson correlation between 
the adjacency matrices of the climate networks in the 
850hPa isobar and in the 500hPa as a function of time. 
Similar to Sec. lIIIBj the upper curve represents the values 
of py for the original networks, the lower curve is after 
the removal of the distance effect, and the middle curve is 
after the removal of both distance and noise effects. As is 



clearly seen from Fig. 7(a) this correlation is significant, 
with an average value of py w 0.6 and small fluctuation 
during time. In all non-equatorial (for both temperature 
and geopotential height) regions we observe similar high 
py values. A further observation from Fig. |7(a)| is the 
increase of 10-20% in pv values after the removal of noise. 

In contrast, in equatorial regions (see e.g. Fig. |7(b)| we 
generally observe lower and more fiuctuative p^ values. 
Another difference between non-equatorial and equato- 
rial regions is the smaller effect of the removal of noise 
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FIG. 7: The correlation between two network adjacency ma- 
trices at different altitudes, 850hPa isobar and 500hPa isobar 
in equatorial and non-equatorial regions, (a) Zone 1 (non- 
equatorial region), and (b) zone 9 (equatorial region). The 
upper curve in each figure, represents the correlation between 
the original networks without removing effects of distance and 
noise. The lower curve represents for the networks after 
removing the distance effect. The mid curve represents p" for 
the networks after removing both, distance and noise effects. 

Based on physical considerations it is reasonable that a 
pair of networks with a larger altitude distance will have 
a smaller similarity. Indeed such a monotonic decreasing 
relation is seen in Fig. |51 Each point in the curves of 
Fig. IS] is an average over different regions, different time 
snapshots, and different altitudes of the correspondence 
py . The different curves indicates that this monotonic 
decrease in the similarity behavior holds both in equato- 
rial and non-equatorial regions and both for temperature 
and geopotential height networks. 

A. Criterion for significant links 

Underlying our supporting arguments for the stability 
of the climate network, there is an assumption. We rely 
on the existence of a sharp boundary between the proper- 
ties of links that result due to noise, L^r, and links that 
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FIG. 8: Average correlations between pairs of networks with 
different altitude distances. The four curves describe two net- 
works which are located in non-equatorial regions (zones 1-6) 
and based on geopotential height or based on temperature and 
two networks which are located in equatorial regions (zones 7- 
9) and based on geopotential height or based on temperature 
measurements. 



result due to real physical dependence, Lp. The set of 
all links is L = Lat U Lp. In this section we will show 
that such a boundary indeed exists, with respect to two 
link properties: (a) the average over all years of the link 
strength , Wi,r and (b) the variability over time of the 
time delays of links STD{Ti^r)- We will later show that 
using both quantities in order to identify the set of real 
links Lp converges to almost the same set of links. 

Our anchor for comparison between the derived Lp and 
Ljv is the distributions of (a) Wi^r and (b) STD{Ti^r) 
for networks based on shuffled data. The shuffling scheme 
is aimed at preserving all the statistical quantities of the 
data, such as the distribution of values, and their auto- 
correlation properties, but omitting the physical depen- 
dence between different nodes (different geographical lo- 
cations). The network properties in such a case are only 
due to the statistical quantities and therefore are similar 
in their properties to false links. To achieve this shuf- 
fling goal, we choose for each node a random sequence 
of y in (the order of the days, d is preserved) fi^. 
Thereafter, the entire construction of the network, based 
on correlations of the shuffled records, is performed. The 
adjacency matrix of the network based on shuffled data 
is denoted as wi_r- The time delay matrix of the network 
based on shuffled data is denoted as ti r- 

(a) Average link strength. In Fig. IH] we compare 
the probability density function (PDF) of Wi_r and wi^r- 
As clearly seen from these figures, the range of possi- 
ble wi,r is extended only over a limited range of values, 
w/,r € [3,4]. Higher values that exist in the PDF for 
Wi^r are missing from the PDF of the shuffled data, and 
therefore are not likely to occur by chance. The cumu- 
lative distribution function (CDF) of wi^r (see insets of 
Fig. [S]) can be regarded as an estimate for the likelihood 
of a Wi,r value to arise by real physical dependence. The 
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98% likelihood level is shaded in the inset of Fig. 9(a)[ 
having Wi,r > 4 for ofF-equatorial regions, and in the 
inset of Fig. [9(b)| Wi^r > 3.6 for equatorial regions. 
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FIG. 10; The distribution of STD{Ti,r) and STD{ti,r) in 
equatorial and non-equatorial regions for networks base on 
temperature measurements at 850hPa isobar, (a) Zone 1 
(non-equatorial region), and (b) zone 9 (equatorial region). 



FIG. 9: The distribution of Wi^r and liJ;,^ in equatorial and 
non-equatorial regions for networks based on temperature 
measurements at 850hPa isobar, (a) Zone 1 (non-equatorial 
region), and (b) zone 9 (equatorial region). 

(b) Variation of the time delay. High variabil- 
ity during different years of the time delays of links 
STD(Tj_r), is also a signature of artificial (random) be- 
havior [4^ . Therefore it serves as another good separator 
between Lp and Lj^. In Fig. [10] we show the probabil- 
ity density function (PDF) of STD{Ti^r) and STD{ti^r)- 
As clearly seen from these figures, the range of possi- 
ble STD{ti,r) is extended over a limited range of values, 
STD{ti^r) e [75,150]. Lower values that exist in the PDF 
of STD{Ti^r) are missing from the PDF of the shufiled 
data, and therefore are not likely to arise by chance. The 
cumulative distribution function (CDF) of STD{ti^r) (see 
inset of Fig. [T0|) can be regarded as an estimate for the 
likelihood of a STD{Ti^r) value to occur by real physical 
dependence. The 98% likelihood level is shaded in the 
inset of Fig.[ini having STD{Ti^r) < 75 for off-equatorial 
regions, and STD{Ti,r) < 80 for equatorial regions. 

Both Wi^r and STD{Ti r) can be used for determining 
a boundary between Lp and Ljv. Convergence to sim- 
ilar Lp and L^r in both criterions can be considered as 



a confirmation that either of these criterions indeed ef- 
ficiently distinct between links that emerge merely due 
to noise and real links. In Fig. [TT] we show a two di- 
mensional PDF of Wi^r and STD{Ti^r)- A large fraction 
of the links evidently have both low values of Wi^r and 
high values of STD{Ti r), which is a typical behavior of 
links that emerge from random behavior. This set of 
links is realized as a sharp local maximum of the PDF 
in the region Wi.r £ [3, 4.5] , STD{Ti^r) e [75, 150] (See 
Figs. Eland Uni). Within this region, Wi.r and STD{Ti^r) 
are not correlated, since the fluctuations are random in 
both axes. Outside this region the mutual local maxi- 
mum of the PDF in both axis are correlated, i.e. larger 
values of Wi^r are paired with lower values of STDiTi^r)- 
The crossover between the two regimes occurs around 
(^i^r,STD{Ti^r)) = (4.5,75). This qualitative behav- 
ior is consistent at all regions (1-9), but the crossover 
point is a bit different for non-equatorial regions (around 
(4.0,75)). 

A further indication that the boundary between Lp and 
ijv is within the region Wi^r G [3.5, 4.5] is the increased 
sensitivity of the stability measure p (r; 9) to the removal 
of noise within this region. Here we indicate explicitly the 
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FIG. 11: (color online). The 2D histogram of links time delay 
variation, STD{Ti,r) and links average strength, Wi^r- This 
2D histogram is for network located at zone 1 and based on 
temperature measurements at 850hPa isobar. 



value of the threshold for noise removal by the second ar- 
gument 9. In Fig. [12] we show the differential of p(r; 9), 
averaged over all values of r, Sp = p{9 + S9) — pi9), 
where 69 = 0.5. We find a sharp local maximum in Sp, 
around Wi.r € [3.5, 4.5]. This sharp maximum is consis- 
tent with 9 crossing the boundary between Lpf and Lp, 
where many of the links related to noise drop off the net- 
work, and causes the average stability p (9) to abruptly 
rise. Such behavior of dp is consistent both in equato- 
rial regions (stars) and non-equatorial regions (circles). 
The response of the sensitivity 6p to further removal of 
links (which mainly belong to Lp) is thereafter reduced. 
In fact, removal of physical links might even result in 
a reduction of the stability (e.g. Sp < 0), as is indeed 
observed for large 9 in the equatorial curve in Fig. [T^l 
In conclusion, a sharp boundary between Ljv and Lp is 
almost certainly identified around Wi.r ~ 4.0 ± 0.5 in 
the networks calculated for all types of data (tempera- 
ture and geopotential height in various altitudes covering 
the troposphere), both in equatorial and non-equatorial 
regions. 



V. SUMMARY 

In summary, we have established the stability of the 
network of connections between the dynamics of climate 
variables (e.g. temperatures and geopotential heights) in 
different geographical regions. The strength of the physi- 
cal connection, Wi^r, that each link in this network repre- 
sents, changes only between 5% to 30% over time. A clear 
boundary between links that represent real physical de- 
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FIG. 12: (color online). The sensitivity, Sp, of the stability 
function p as a function of the threshold 9. 



pendence and links that emerge due to noise is shown to 
exist. The distinction is based on the high link strength 
Wi^r and on the variability of time delays STD {Ti,r). 

Beside the stability of single links, also the hierarchy 
between the link strengths is preserved to a large extent. 
We have shown that this hierarchy is partially due to the 
two dimensional space in which the network is embedded, 
and partially due to physical coupling processes. More- 
over the contribution of each of these effects, and the 
level of noise was explicitly estimated. The spatial effect 
is typically around 50% of the observed stability, and the 
noise reduces the stability value by typically 5%-10%. 

The network structure was further shown to be con- 
sistent across different altitudes, and a monotonic rela- 
tion between the altitude distance and the correspon- 
dence between the network structures is shown to exist. 
This yields another indication that the observed network 
structure represents effects of physical coupling. 

The stability of the network and the contributions of 
different effects were summarized in specific relation to 
different geographical areas, and a clear distinction be- 
tween equatorial and non-equatorial areas was observed. 
Generally, the network structure of equatorial regions is 
less stable and more fluctuative. 

The stability and consistence of the network struc- 
ture during time and across different altitudes stands in 
contrast to the known unstable variability of the daily 
anomalies of climate variables. This contrast indicates 
an analogy between the behavior of nodes in the climate 
network and the behavior of coupled chaotic oscillators 

A future outreach of our work can be a mapping be- 
tween network features (such as network motifs) and 
known physical processes. Such a mapping was previ- 
ously shown to exist (29j between an autonomous clus- 
ter in the climate network and El-Niiio. Further struc- 
tures without such a climate interpretation might point 
towards physical coupling processes which were not ob- 
served earlier. 
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